%% reconstruct_from_weights.m
%
% This function estimates the source space activity from the SPM object and
% the weights of the beamformer. It is to be used when the Y_evoked data
% has not been saved in order to conserve disk space.
%
% dat_source = reconstruct_from_weights(source_recon_results,D,voxels2use,trials2use)
%
% Inputs:
% source_recon_results : the output object from osl_beamformer_fuse_scalar
% D : the input SPM object.
% voxels2use : the index of the desired voxels (default is all).
% trials2use : the index of the desired trials (default is all).
%
% Use this with caution as it has not been exhaustively debugged.
%
% HL 07.09.12

function dat_source = reconstruct_from_weights(source_recon_results,D,voxels2use,trials2use)

%Set-Up
if nargin < 4, trials2use = 1:length(source_recon_results.trialtypes); end;
if nargin < 3, voxels2use = 1:size(source_recon_results.grid.pos,1); end;
if nargin < 2, D = spm_eeg_load(source_recon_results.S.D) ; end;

Nsamples = length(source_recon_results.times);
Ntrials  = length(trials2use);   
Nvoxels  = length(voxels2use);

% Apply bandpass filter to good channels
allchans=[];
for ff =1:length(source_recon_results.S.modalities)
allchans=sort([allchans find(strcmp(D.chantype,source_recon_results.S.modalities{ff}))]);
end
dat_sens = D(setdiff(allchans,D.badchannels),:,:);
for i=1:size(dat_sens,1)
    for j=1:size(dat_sens,3)
        dat_sens(i,:,j) = bandpass(dat_sens(i,:,j) ,source_recon_results.freq_range,D.fsample);
    end
end

% Apply Normalisation to different channel types.
for ff = 1:length(source_recon_results.S.modalities)
    sensind=find(strcmp(D.chantype(setdiff(allchans,D.badchannels)),source_recon_results.S.modalities{ff}));
    dat_sens(sensind,:,:) =  dat_sens(sensind,:,:)/source_recon_results.normalisation{ff};
end

% Estimate activity at each voxel
dat_source = zeros(Nvoxels,Nsamples,Ntrials);
for vox=1:Nvoxels
    for tri=1:Ntrials
        dat_source(vox,:,tri) = source_recon_results.weights(voxels2use(vox),:)*dat_sens(:,:,trials2use(tri));
    end
end